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Abstract 

Schrodinger equation with given, a priori known current is formulated. A non-zero current 
density is maintained in the quantum system via a subsidiary condition imposed by vector, local 
Lagrange multiplier. Constrained minimization of the total energy on the manifold of an arbitrary 
current density topology results into a non-linear self-consistent Schrodinger equation. The applica- 
tions to electronic transport in two-terminal molecular devices are developed and new macroscopic 
definition of a molecular current-voltage characteristic is proposed. The Landauer formula for the 
conductance of an ideal one-dimensional lead is obtained within the approach. The method is 
examined by modeling of current carrying states of one-dimensional harmonic oscillator. 
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I. INTRODUCTION 



The remarkable miniaturization of conventional microelectronic devices has been wit- 
nessed over the last decades. If the miniaturization trend is to continue, elements of micro- 
electronic circuits, e.g. transistors and contacts, will soon shrink to single molecules or even 
atoms. In the past few years, there have been considerable advances towards fabrications 
and experimental studies of mono-molecular electronic devices [jl|. This activity has been 
largely spurred on by development of experimental techniques to form atomic scale electri- 
cal contacts such as scanning tunneling microscopy M and mechanically controllable break 
junctions [0]. Measurements of current- voltage (I-V) characteristics have been recently per- 
formed on a single benzene- 1,4-dithiolate molecule ||, buckminsterfullerene ||, individual 
atoms H, nanoclusters Q] and carbon nanotubes ||. These initial observations of con- 
duction in molecules were followed quickly by experimental demonstration of the simplest 
devices: molecular diodes f9|, switches [jlOj and memory elements |11| although molecular 
equivalent of a transistor is yet to be discovered. 

Several theoretical approaches have been recently developed to describe electronic con- 
ductivity of molecules [0. One of the widely used methods is based on the combination of 
extended Hiickel model or density functional theory to treat the molecular electronic struc- 
ture, an electronic reservoir model for the contacts and linear response approximation to 
compute the electronic current as a function of the applied voltage [|13|, [b|, [15], [16], [T7). The 
most recent applications used Lippmann-Schwinger or non-equilibrium Green's functions 
formalisms plus "jellium" model for the contacts and for the molecular-contact interactions 



fl8| , pT9 1 . In all these approaches, the system is partitioned into three parts: molecular wire 
and two electron reservoirs to represent metallic contacts. The use of the model dependent 
reservoir-molecule interaction as one of the most pivotal ingredient diminishes the ability of 
these methods and prevents from a development of predictive and quantitative techniques 
for the electronic transport problem. Recent density functional theory based calculation 
demonstrates that a molecular wire is bonded to a gold contact by the thiol bridge which is 



known to be the strong covalent interaction [2D]. The empirically treated coupling between 
metal contacts and the molecular wire, which serves as the interaction for the linear response 
calculations and represents the kernel of the corresponding Lippmann-Schwinger equation, 
is at least strong as intermolecular interaction. Recent papers also indicate the exact nature 



2 



of the molecule-contact bonding is critical in predicting the correct order of magnitude for 
the current with applied bias 

The miniaturization shrinks the size of the system to a level when the problem can be 
addressed by ab initio electronic structure methods. As an example consider a standard 
prototype molecular device system with a single benzene- 1,4-dithiolate molecule covalently 
bonded to two gold electrodes M. The problem can be reduced to the system where the 
molecule is bonded to two small Au clusters [|17] . A gold cluster with reasonably small edge 



influence on the molecular-surface bonding contains 10 atoms [20]. With a single benzene- 
1,4-dithiolate molecule as the molecular wire the contact- wire-contact system consists of just 
32 atoms. Systems of such size can be routinely simulated by standard quantum chemistry 
codes. The main obstacle for the application of the first principle electronic structure meth- 
ods is not the size of the system but modeling of the electronic current and interpretation 
of the applied voltage. Suppose that we have solved Hartree-Fock or Kohn-Sham equations 
for a molecular device. The ground state wave function is real in absence of a magnetic field 
and therefore the current density computed from this wave function is zero everywhere in 
the sample. Another important issue is an interpretation of the applied voltage. Virtually 
all present theoretical approaches are pivotally relied on simplified, noninteracting electrons 
type models for the metal contacts where the external voltage bias can be straightforwardly 
associated with the Fermi energies shift between left and right electronic baths. One of 
the reasons for use electronic reservoir models to represent the metallic contacts is principle 
theoretical difficulties related to a rigorous quantum mechanical definition of the applied 
voltage. 

In this paper, we advocate an alternative approach where the use of the electronic reservoir 
to model contacts is completely avoided from the very beginning and a quantum system with 
current is treated effectively as a closed system. We begin by describing the variational 
principle for a quantum system with current. We then derive the Schrodinger equation 
with probability flux which provides exact wavefunctions on the manifold of a given current 
density. We next specify two-terminal molecular device steady current subsidiary condition 
and develop the Schrodinger equation for a two-terminal molecular device. We discuss the 
use of the energy ballance condition to compute current- volt age characteristics of quantum 
systems. We then derive the Landauer formula for non-interacting charge carriers within 
the formalism. Finally, we demonstrate the salient features of the approach through the 
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numerical solution of the self-consistent Schrodinger equation for one-dimensional harmonic 
oscillator with constrained current. 

II. SCHRODINGER EQUATION WITH PROBABILITY CURRENT 
A. Variational principle 

We start with the nonrelativistic Hamiltonian for a particle in the external scalar potential 
v(r) (atomic units are used throughout the paper unless otherwise mentioned): 

H = -\v 2 + v{v). (1) 

We define density and paramagnetic current density as: 

p(r) = V(r)V*(r) (2) 

j(r) = Im{V>*(r)W(r)} (3) 

where ip(r) is not generally an eigenvector of the Hamiltonian H . We will specify an eigen- 
problem for the ip{r) in the next section. The physical current is equivalent to the param- 
agnetic current densities in absence of an external magnetic field. 

An eigenvector for a bound state of a real Hamiltonian can be always made real by 
a change of the phase of the wave function. The Hamiltonian H is real therefore if we 
just solve standard Schrodinger equation with the Hamiltonian H, we will obtain real wave 
function ip(r) and the current density j(r) will be automatically zero. The standard way 
to get a complex wavefunction out of a real Hamiltonian is to assume certain boundary 
conditions. These boundary conditions could be an assumption of an incoming plane wave 
in scattering theory or periodic boundary conditions for a translationally invariant system. 
To solve a scattering problem for a molecular device is possible only at the cost of major 
unphysical approximations, e.g. via considering the molecular device as a singe scattering 
center. Another approach is to assume periodic boundary conditions. The periodic boundary 
conditions prevent one from studying any situations where the change of chemical potential 
across the system is of finite magnitude because the chemical potential must be also periodic 
in space. The local chemical potential is the decreasing function in the direction of the 
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current and therefore the system can not be treated as translationally invariant at least 
along the current axis. 

We have developed an alternative approach where a current carrying quantum system is 
treated effectively as a closed system with a subsidiary condition for the generally discontin- 
uous current density distribution. We start our derivation with the assumption that current 
carrying states of an open quantum system can be described in terms of wave functions of an 
effectively closed quantum system if the current is maintained via explicit constraint. This 
assumption is a simplified, particular case of the standard assumption in statistical mechan- 
ics: a system coupled to a thermal bath can be described with statistical ensembles of a 
closed system in which the system-bath interaction is not explicitely considered, but enters 
only via Lagrange multipliers and subsidiary constraints. It has been recently demonstrated 
that a steady current open mesoscopic system can be described in terms of an ensemble 
of carrier states in the closed mesoscopic device itself if the steady current is additionally 



constrained [^TJ p2j] . 

We require that the current density is constrained to be specified function I(r). It results 
in the subsidiary constraint equation: 

j(r)-I(r)=0. (4) 

The expression (|j) generally can not be obtained by differentiating of a functional of ip and ip* 
with respect to r. It means that the constraint eq.(^j) is nonholonomic [23|. The variational 



problem is to minimize the total energy with respect to ip subject to the nonholonomic 
constraint for the current density @). Two types of treatment are possible: we can either 
choose in which the constraint is implicit or we can impose the constraint explicitely 
by using the Lagrange multiplier. We will follow the second approach which was recently 
suggested by Kosov and Greer |24[] , elimination the part of the wavefunction by imposing the 



explicit constraint. Although the constraint formulated as eq.([|) is vector and nonholonomic 
it can be included into a variational functional via pointwise, vector Lagrange multipliers 
121. 



a r 



Within the Lagrange multiplier approach the constraint is explicitly achieved via intro- 
duction of the auxiliary functional: 

= (ij>\H\ij>) - E «V#> - 1) - J a(r) ■ (j(r) - I(r)) dv . (5) 
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The first two terms are standard in the variational derivation of the Schrodinger equation 
with the first giving the total energy and the second is introduced to maintain the normal- 
ization of the wave function. The third term with the vector Lagrange multipliers a(r) has 
been introduced to impose the subsidiary constraint for the current density. 

A use of the variational principle can be justified only for closed quantum systems. The 
choice of I(r) can not be obviously completely arbitrary to maintain the closed system 
boundary condition. The localized in space quantum system will be considered closed if, 
first, there are no in- and out- probability flows in the system; and, second, if the wave 
function ip is an eigenvector of a Hermitian operator (it prevents the complex eigenenergies 
and therefore effectively closes the system). We will demonstrate in the next section that 
both these conditions are satisfied within our approach. 



B. Self-consistent Schrodinger equation 

The minimization of the total energy subject to the subsidiary constraint for the current 
density (f|), i.e the variation of the auxiliary functional £l[if>], results into the following non- 
linear Schrodinger equation: 

|a + v(r) - 1 [a(r), V] + ) ^(r) = E^{v) . (6) 

The anti-commutator term, 

i[V,a(r)] + = l(Va(r)+2a(r)V), (7) 

is the additional completely imaginary potential arising directly from the constraint on the 
current density (consult with appendix A for the corresponding functional derivatives ). This 
additional constraint potential forces the wave function ^(r) to be irremovably complex and 
insures that ^(r) produces the required current density I(r). 

A physical interpretation of the Lagrange multiplier a(r) can be made more explicit if we 
re- write eq.(^J) in the following form: 



I(_, V -a(r)) 2 + ia(r)-a(r) + W (r) 



V(r) = E^{v) . (8) 



It can be now established that the current constraint has introduced the terms which are 
mathematically equivalent to an external vector potential. But unlike an external vector 



potential, the Lagrange multiplier a(r) does depend upon the wave function ip(r) and eq.@ 
becomes nonlinear and must be solved self-consistently. The eigenvalue problem eq.(|8]) is 
generated by Hermitian operator since the Lagrange multiplier a(r) is real. Being a solution 
of the Hermitian eigenproblem the eigenenergy E is real and the wavefunctions ip form a 
complete set in Hilbert space. Within our approach we deal with the localized in space 
system with the wave functions generated by the hermitian eigenproblem (|8]). Such system 
is by definition closed and therefore the use of the variational principle is formally justified. 

The Schrodinger equation with current (||]) is not yet in a form allowing for a solution to 
be found as the Lagrange multiplier a(r) is not known yet. We next seek a solution for the 
Lagrangian multiplier a(r). In the developing of this method we have found out that the 
most significant source of numerical instabilities in the solution of the nonlinear Schrodinger 
equation (|]) arises from the determination of the Lagrange multiplier a(r). The described 
below procedure obtaining of the Lagrange multiplier has been found to be the most robust. 
Multiplying eq. (||) from the left on the ip( T )* an d subtracting from the obtained equation 
its own Hermitian conjugate, yields 

V-j(r) = V-(a(r)p(r)). (9) 

Equation (^) can be resolved for a(r) with the assumption that both the density p(r) and 
the current density j(r) decay to zero at infinity. In this case the direct integration of eq.@ 
results into the following expression for the Lagrange multiplier a(r): 

a(r) = || (10) 

Substitution of the constraint for the current density eq.(|j) into eq.(|T0D yields the simple 
expression in which the Lagrange multiplier depends upon the wave function only via the 
density in the denominator: 

a(r) = ^ (11) 



There is a direct similarity between equation for the Lagrange multiplier (|T0|) and the 
definition of the velocity of a quantum trajectory within de Broglie-Bohm causal description 



of quantum mechanics |j25 |. Within the de Broglie-Bohm dynamics the velocity u(r, t) of a 



trajectory at a given point is computed as f26fl: 



u < r ^M- (i2) 



therefore, at least for stationary systems the vector Lagrange multiplier a(r) can be directly 
related to the velocity for a corresponding de Broglie-Bohm trajectory. By this analogy we 
can make the physical interpretation of the Lagrange multiplier more explicit. In the de 
Broglie-Bohm quantum dynamics to maintain the constant flux at the given density profile 
one needs to adjust the local trajectory velocities. The vector Lagrange multiplier appears 
to play the simular role in the Schrodinger equation with constrained current. Likewise for 
a point with the non-zero trajectory velocity, the current density has finite value at a point 
with the non-zero vector Lagrange multiplier. And likewise the de Broglie-Bohm velocity, 
a product of the vector Lagrange multiplier on the probability density yields the current 
density. 

With this choice of the Lagrange multiplier ([TT|) we obtain the Schrodinger equation with 
the current I(r) in the self-consistent and closed form: 

1 / I(r)\ 2 l/ 2 (r) 



. . , , , _ V(r) = i#(r). (13) 

2 V p{r)J 2p 2 (r) 

The derivation of the self-consistent Schrodinger equation eq.(|T3"D composes one of the central 
results of the paper. 



III. SCHRODINGER EQUATION FOR A TWO TERMINAL MOLECULAR DE- 
VICE 

The desire to develop the basic Schrodinger equation for calculations of I — V character- 
istics of two terminal molecular devices is the impetus for the present study. To achieve this 
aim, we specify the molecular wire constraint on the current density: 

J dydzj x {r)=I x 6(x-L,R-x) , (14) 

where I x is the steady current through the molecular device. The step function, 

6(x- L,R-x) = 6(x- L)6(R-x) , (15) 

equals to 1 if L < x < R and zero otherwise. L and R are the left and right boundaries 
for the molecular device. Within this description, net current flow is aligned along the 
and there are no in- and out flow: the current starts at the left boundary L and 
completely absorbed at the right R. With this constraint it is not required to fix the full 
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vector j(r) everywhere in the sample. Only the net current flow across a cross section 
J dydzj x (r) is constrained, and this quantity is readily available experimentally. This is a 
simple geometric arrangement to specify current flow for the two terminal molecular device 
and can be extended to more complex device topologies, e.g. for the three terminal devices 
where three tips are connected to a single molecule JZ7 . 

Likewise, as for the derivation of the exact Schrodinger equation in the preceding section, 
with this choice of the constraint ( |14|) the minimization of the auxiliary functional @ now 
yields the Schrodinger equation for the steady current two terminal molecular device: 



1 A / N 1 

-2 A + V ^-2l 



a x (x), 



d_ 

dx 



if)(r) = Eip{v) 



(16) 



The expression for the Lagrange multiplier given by eq. ([11]) can be straightforwardly 
extended to the steady current molecular device, with the result is 

I, 



a x (x) 



Pyz\ x 



-6{x-L,R 



x 



(17) 



where we have introduced the quantity 



Pyz{x) 



dydzp(r). 



(18) 



With the expression for the Lagrange multiplier a x (x) in hands, we can complete the 
Schrodinger equation for the two-terminal molecular device with steady current I x {x) 

(x — L, R — x) 



1a 

2 



d_ 

dx 



(19) 



Some aspects of eq.([TJ) deserve special discussion. The imaginary potential which enforces 
the wave function to produce steady current I x depends upon the x coordinate only. The 
range of this potential is restricted by the molecular device boundaries. One of the terms in 
this potential is proportional to the derivative of the step function, 



^-9{x -L,R-x) = 5(x-L)- S(x - R) 
ox 



(20) 



and therefore yields the singular imaginary 5-function potential on the device boundaries. 
It is known from the standard quantum mechanics textbooks p8f that a 5-function poten- 
tial results in the discontinuity of the first derivative of the corresponding solutions of the 
Schrodinger equation. Likewise, the imaginary component of the eigenfunction ^>(r) of the 
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Schrodinger equation for a molecular device (|T9|) has the discontinuous first derivative at 
the points x = L and x = R. 

The final comment regarding the computation of the I-V characteristic of a two-terminal 
molecular device is in due order. The current carrying states obtained as the solutions of 
the molecular device Schrodinger equation (|T^) must now be associated with current-voltage 
characteristics available experimentally. The rigorous definition of the applied voltage is a 
controversial issue if one does not invoke to noninteracting electron reservoirs to represent 
contacts. We propose the following macroscopic definition of the applied voltage which has 
the direct relation with the quantities measured experimentally. The energy increase from 
the applied voltage U is computed by the numerical integration with an assumption of linear 
voltage drop between the left and right boundaries: 

R 

AE d (U) = J dx(R - x)p yz (x) . (21) 

L 

We compute the energy of the non-active, i.e. with zero current, molecular device and 
compute the total energy of the same molecular device at an experimentally given current 
I x . Then the energy increase for the establishment of the current carrying state, i.e. the 
energy difference 

AE d (I x ) = E d (I x ) - E d (I x = 0) (22) 

is calculated. Because of the energy ballance condition the AE d {U) must be equal to the 
AE d (I x ) and it results into the following equation for the applied voltage bias: 

U = ^— ^ AE d (I x ) . (23) 

/ dx(R - x)p yz (x) 

L 



The eigenenergy E of the Schrodinger equation (|19f) should not be confused with the energy 
of molecular device. The eigenenergy E serves only as the Lagrange multiplier to enforce 
orthogonality of the wave function. The energy of the molecular device can be obtained as 
an expectation value of the real physical Hamiltonian H with the integration restricted by 
the molecular device region: 

R 

E d (I x ) = J dydz J dxr(r) (~A + v(r)j ^(r) , (24) 
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where ip(r) is the solution of the Schrodinger equation (|T9|) with the given I x . Since the 
wavefunction ^(r) is complex and the integration (^) is not performed over whole space, 
it is not obvious from the expression fl24| ) that the device energy Ed(I x ) must be real. We 
will show now that the device energy E d (I x ) is real if the wavefunction ip(r) produces the 
constrained current density flUp. To prove this we subtract from eq.(^) its own complex 
conjugate and after a bit of algebra arrive to the following expression: 

R 

E d (I x ) - E d (I x y = -i J dydz J dxW ■ j(r) = -i I dSn ■ j(r) . (25) 

The length of the quantum wire is aligned along the x— axis and the surface integral does 
the integration over the surface of the rectangular box with one y — z face at x = L, another 
y — z face at x = R and any of the x — y or x — z faces are removed at the infinity. Within 
our constraint, there is a positive current flowing into the y — z face at x = L and exactly 
the same positive current flowing out of the y — z face at x — R. There is no current flow 
out of the region in any of the x — y or x — z faces. Therefore the surface integral vanishes 

j) dSn-j(r) = (26) 

and the device energy Ed(I x ) defined via eq.(|2~3]) becomes real. 



IV. EXAMPLE CALCULATIONS 



A. Conductance of an ideal one-dimensional lead 



In order to establish connection with standard approaches and test the method we derive 
the Landauer formula |?9| using the Schrodinger equation for current carrying states eq.(|T3D 
as the starting point. Consider noninteracting electrons moving along a one-dimensional 
wire without any scatterer. The electronic density does not depend upon x yielding the 
following Schrodinger equation with steady current /: 



Id 2 .Id 

2 dx 2 p dx 



if> k (x) = E k ip k (x) . 



(27) 



Equation (|27D has the plane-wave solution 



i'ki.x) ~ exp(ikx) 
11 



(28) 



with eigenenergy 

_ k I k, , . 

E k = -- — . (29) 
z p 



The similar dispersion relation as given by eq.(P9"D has been obtained within the same ap- 
proach for uniform electron gas with applied direct current [p4j] . For the small voltage and 



temperature transmissions of electrons take place only through the states in close vicinity 
of the Fermi level. There are two vectors k, one is positive and one is negative, which 
correspond to the Fermi energy: 

= ? + \/ 2 ^ + * ' (30) 

k - =L p-f^W' <3i) 

Following the Landauer approach we associate the voltage drop (at the zero temperature) 
with the gap between the singe-electron energies of the electrons moving in the direction of 

the current (k > 0) and electrons moving in the opposite direction (k < 0): 

P k 2 

f = f- T . (32) 



Substituting the expressions for the k + and eqs. (|30| , |3TD into the Landauer's definition of 
the applied voltage (|32D, we readily find 

/ = GU , (33) 

with the conductance 

G = 9 . (34) 

2 ] /2E F + (I)' 

Then, the use of the standard relation between the Fermi momentum kp and the electronic 
density p = k F /ir along with the small current expansion of the square root eq.(^) results 
into the Landauer formula for a single transport channel of an ideal one- dimensional lead 

[3g : 



G = JL. (35) 

If we take into account the spin degeneracy of real electrons, the conductance becomes 
multiplied by 2. By accounting the degeneracy and converting the conductance eq.(|35|) to SI 
units the standard value G = e 2 /(irh) = (12.9 KQ)" 1 is obtained. Re-deriving the standard 
result we demonstrate that in the low current and zero temperature limits our method is 
equivalent to the traditional approaches which are based on the Landauer formula. 
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B. Numerical solution of the Schrodinger equation with current 



We demonstrate applications of the method with numerical examples. The numerical 
solution of a self-consistent Schrodinger equation is often performed within a finite basis set 
expansion for the wave function with subsequent diagonalization of the Hamiltonian matrix 
until the self consistency is reached. We have found out that it is computational convenient to 
form the complete, real basis set from the solutions of the zero-current Schrodinger equation: 

f_lA + V (r)V(r) = e M ^(r). (36) 
Then we expand the wave function of a current currying state on this basis 

^(r) = ^C M ^(r), (37) 

the coefficients C M must be complex for the wave function which produces a non-zero current. 
With the expansion (|37|) we arrive to the following complex, Hermitian eigenproblem: 

£ (e^ - mi>„ m (I)) C, = EC V , (38) 

where the current dependent matrix element w VIM (T) is computed in appendix B and is given 
by the following integral: 

«V(I) = i Jdv^r- (0 m (r)V0;(r) - ? ;(r)V0 m (r)) . (39) 

The matrix elements w Ufl (i) are real and depend on the expansion coefficients via the 
density p(r) thereby making the eigenvalue problem ( |38"D nonlinear. The eigenproblem ( |38| ) 
is Hermitian because the matrix element w VfM (l) changes its sign when we swap \i and v 
indeces: 

w vlt (J) = -iiv(I) . (40) 

The numerical calculations were carried out for the one- dimensional system. With an 
eye on molecular device simulations we solve the self-consistent Schrodinger equation for the 
step function constraint as a primary example: 

j(x) = I 0(x-L,R-x) (41) 

The external potential was taken in the harmonic oscillator form and we put L = —2.0 
and R = 2.0 as left and right boundaries for the current carrying part of the system in all 
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our calculations. Our basis set was built up from the 25 eigenvectors of the zero current 
Hamiltonian. The self-consistency was assumed to be achieved when the root mean square 
deviation of the wave function vectors of the two subsequent iterations was less then 10~ 8 . 

First, we plot the Lagrange multiplier as a function of the coordinate on the figure 1. The 
Lagrange multiplier is discontinuous at the system boundaries. The Lagrange multiplier is 
vanishing out of the system region and it is the inverse function of the density p(x) inside 
the system. As it can be expected from the analogy with the de Broglie-Bohm quantum 
hydrodynamics (|12|) (the lower the density, the larger the local velocity) the Lagrange mul- 
tiplier reaches its maximum values in the regions of depleted density and has the minimal 
value at the point of maximal density. 

In our next example we studied the deviation of the current carrying wave function ip(r) 
(I = 0.1) from the zero-current ground state wave function </>o( r )- The real component is 
depicted on the upper panel of the fig. 2. The shape and magnitude of the real component 
do not undergo significant transformation with the establishment of the current carrying 
state. Although the admixture of the higher terms in the expansion (^) results in the small 
damping oscillations on the tails of the real component of the wave function. The imaginary 
part of the wave function ( plotted on the lower panel of the fig. 2) is identically zero when 
the current is vanished. The imaginary component reaches its negative minimum at the left 
boundary L = —2.0 of the current carrying part of the system and then it discontinuously 
changes the sign of the first derivative and undergoes almost linear increase until it reaches 
its maximum at the right boundary R = 2.0. When the coordinate x goes out of the non-zero 
current region the imaginary part of the wave function rapidly decays to zero value. 

The establishment of the imaginary part of the wave function directly indicates the devel- 
opment of the current carrying state. We studied transformations of the imaginary compo- 
nent as a function of the increased current. The imaginary components of the wave functions 
are plotted on the fig. 3 for three representative steady currents: / = 0.01, / = 0.5 and 
/ = 0.1. At the very small current I = 0.01 the imaginary wave function is almost zero 
everywhere except the only small negative minimum at the left boundary and with the posi- 
tive maximum at the right. The magnitudes of the imaginary component at the extremums 
are nearly linearly proportional to the value of the constrained current. The slope of the 
imaginary component, i.e. the derivative, is always positive and almost constant between 
the left and right boundaries of the system. The imaginary part becomes zero when the 
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real component reaches its maximum, i.e. at the point x = 0.0. Although the total wave 
function is real at the point x = 0.0 the current is not vanished because the derivative of 
the wave function remains complex. 

Now we turn to the calculations of a "device" energy required to establish the current 
carrying state. Being the Lagrangian multiplier to maintain the wavefunction orthonormal- 
ity, the eigenenergy E of the eigenproblem ( |3"ED is not directly related to the energy of the 
system. The energy increase due to the establishment of the current carrying state is defined 
via the space restricted expectation value of the Hamiltonian (|24]) : 

R R 

E d (I) = J dxr(x)(-\^ + v(x)^(x) = Y,C;C u e v j dx <j>p(x)Mx) ■ (42) 

L ^ v L 

The table 1 contains the energy E d {I) and the energy difference between the nonactive, i.e. 
zero-current, and current I carrying state as a function of the increased current. The "device" 
energy is gradually increased with increasing current. As it can be deduced from the table 1, 
we need to pump into the system ~ 0.14 a.u. of energy to establish the current / = 0.1 a.u. 
carrying state. A simple fit of the data from the table 1 provides only the quadratic and 
cubic dependences of the "device" energy upon the steady current (for I < 0.1): 

E d (I) = E d (I = 0) + 17.4/ 2 - 36.5J 3 , (43) 

thereby demonstrating that the term linear in the current has the minor importance for the 
energetics of the current carrying quantum system and the main energy contribution comes 
from the I 2 term. This energy dependence is generally agreed with the current-density 
functional theory of inhomogeneous interacting electron gas where the current dependent 
correction to the standard exchange- correlation energy functional is also mainly proportional 
to I 2 flp. 



V. CONCLUSIONS 



In this paper, we have given a variational formulation of the Schrodinger equation with 
non-zero current. The Schrodinger equation with current is derived via a constrained min- 
imization of the total energy with a subsidiary condition for the current density. The sub- 
sidiary condition for the current density is maintained during the course of variation by a 
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vector, pointwise Lagrange multiplier. An explicit elimination of this Lagrange multiplier 
results into the closed, self-consistent form of Schrodinger equation with current. We showed 
that the current carrying states are eigenvectors of the complex, Hermitian operator. The 
formulation has been developed for general current density topologies and then specified for 
the case of a two-terminal molecular device. We showed how the energy ballance condition 
can be used for rigorous computation of I-V characteristics of molecular devices: the energy 
increase due to the establishment of a current carrying state is associated with the applied 
voltage bias. We demonstrate the salient features of the approach by re-deriving the Lan- 
dauer formula for the conductance of an ideal one-dimensional lead and through a solution 
of the one-dimensional Schrodinger equation with fixed current density. The complex, Her- 
mitian matrix was formed on a basis of the eigenvectors of the zero-current Hamiltonian. 
The numerical solution was achieved via the self-consistent solution of the complex, Hermi- 
tian eigenproblem. We found out that a basis set formed from the zero current Hamiltonian 
eigenvectors provides a robust algorithm for the self-consistent convergence. Our method 
has the great compatability with standard electronic structure methods which are also based 
on a variational principle e.g. Hartree-Fock or configuration interaction approximations. 
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APPENDIX A: FUNCTIONAL DERIVATIVES OF THE CONSTRAINT FUNC- 
TIONAL 



In this Appendix we compute the functional derivative of the constraint functional. The 
constraint functional has the following form 

AM = J dr a(r) • (j(r) - I(r)) . (Al) 

The direct variation of the constraint functional with respect to the ip*(r) yields: 



Si/j*(r) J 5j(r') #*(r) J v ' 5ip*(r) 

= I J rfr'a(rO-^y{^(r')V^(r')-^(r')V^(r')} 

= ia(r) ■ V^(r) - 1^ | dr' a(r') • ^(r')WV) 

= la(r) ■ V^(r) - 1^ | dr' V • {a(r>(r>*(r')} 



= l a (r) • V^(r) - ^ dS n • a(r>(r')^(r') + • {a(r)^(r)} 

= Ia(r) • W(r) + ly • {a(r)^(r)} = I [a(r), V] + ^(r) 



APPENDIX B: MATRIX ELEMENTS OF THE ANTI-COMMUTATOR [a(r), V] + 

The computation of the matrix element of the anti-commutator is presented below. The 
notation V a means that the gradient acts on the Lagrange multiplier a(r) only. 



(<j) n \ [a(r), V] + |0 m ) = (0n| V a • a(r) + 2a(r) • V|0 m ) 

= J dv 0;(r)0 m (r)V • a(r) + 2 J dr a(r) • <ft(r) V0 m (r) 

= y drV-{a(r)0;(r)0 m (r)}- | dr a(r) • V{C(r)0 m (r)} + 2 | dr a(r) • C(r)V0 m (r) 
= jdSn- a(r)0;(r)0 m (r) - | dr a(r) • (0 m (r)V0;(r) - #(r)V0 m (r)) 
= - | dr a(r) • (0 m (r)V0;(r) - 0;(r)V0 m (r)) 
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Table 1. The energy of the active, i.e current carrying L < x < R, part of the system as a 
function of the steady current. AE d = E d (I) — E d (I — 0) is the energy increase to establish 
of the current / carrying state. 



/ 


E d {I) 


AE d (I) 


0.00 


0.4977 


0.0000 


0.01 


0.4994 


0.0017 


0.02 


0.5043 


0.0066 


0.03 


0.5124 


0.0147 


0.04 


0.5233 


0.0256 


0.05 


0.5368 


0.0391 


0.06 


0.5526 


0.0549 


0.07 


0.5706 


0.0729 


0.08 


0.5905 


0.0928 


0.09 


0.6122 


0.1145 


0.1 


0.6356 


0.1379 
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Figure 1. The Lagrange multiplier a(x) is plotted as a function of the coordinate x. 
Calculations were performed for the step function current density distribution j = 0.1 9(x — 
2,2 -x). 




Figure 2. The spatial dependences of real and imaginary components of the current carrying 
wave function. Calculations were performed for the step function current density distribution 
j = 0.1 9(x — 2, 2 — x). The dashed line on the upper panel is the ground state wave function 
of the system with zero current. 
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Figure 3. The spatial dependences of imaginary components of the wave functions are 
plotted for different values of the current. The current density is constrained to be the step 
function j — I 9(x — 2, 2 — x). The solid line corresponds to / = 0.01, dotted to I = 0.05 
and dashed to / = 0.1. 
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